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Abstract 

A formalism for the numerical integration of one- and two-loop integrals is pre- 
sented. It is based on subtraction terms which remove the soft, collinear and some of 
the ultraviolet divergences from the integrand. The numerical integral is performed 
in the Feynman parameter space, using a complex contour deformation to ensure ro- 
bust convergence even in the presence of physical thresholds. The application of the 
proposed procedure is demonstrated with several one- and two-loop examples. An 
implementation in the program NICODEMOS is publicly available, which currently in- 
corporates only one-loop functionality, but an extension to two-loop cases is planned 
for future versions. 



1 Introduction 



For many production processes at colliders and other particle physics observables, the inclu- 
sion of higher-order radiative corrections is essential to match the experimental precision. 
Much progress towards efficient calculational techniques has been made over the last decades 
by many authors. However, the computation of loop diagrams with many different mass 
scales, many external legs, or more than one loop remains a difficult task. In particular, 
when going beyond the one-loop level, it is known that the loop integrals can in general not 
be solved analytically. Nevertheless, by using specialized semi-numerical techniques, several 
complete two-loop calculations have been carried out, for example for electroweak precision 
observables [Il[2] and QCD corrections to the production of gauge bosons [5], top quarks |1] 
or Higgs bosons [S] at the Large Hadron Collider (LHC). However, the methods employed 
in these papers have been tailored to the problem at hand and cannot be applied easily to 
other situations. 

Fully numerical techniques offer an alternative and potentially more flexible approach to- 
wards complex loop calculations. For a purely numerical method to become viable, two main 
obstacles need to be overcome: extracting the ultraviolet (UV) and infrared and collinear 
(IR) singularities from the integral; and ensuring robust and efficient convergence of the 
numerical integration. Several powerful methods have been proposed: 

• Sector decomposition starts from the Feynman parametrization of a loop diagram. 
It isolates the physical singularities through appropriate mappings of the integration 
region There are additional integrable singularities inside the integration region 
for a diagram with physical cuts. While being formally integrable, such singularities 
cannot be handled by standard numerical integration algorithms. To avoid the problem 
and improve numerical stability one can deform the integration contours in the complex 
plane [7]. 

• The Feynman integrals can be transformed into Mellin-Barnes representations. An 
algorithm for the extraction of UV and IR poles has been developed based the residue 
theorem jSJlS] . Again, the presence of physical cuts leads to bad convergence behavior 
of the numerical integrals, which can be improved by variable mapping and contour 
deformations [TO]. 

• In a third approach the UV and IR singularities are removed from the integrand by 
suitable subtraction terms |lll[12]- These subtraction terms are designed such that 
they can be easily integrated analytically and added back to the final result. As for 
the previous two methods, integrable singularities due to physical thresholds may be 
treated with the help of contour deformations (T2] . 

The method based on subtraction terms has so far been applied only to one-loop calculations. 
This article presents a procedure for extending this idea to two-loop integrals. For this 
purpose the IR subtraction terms are constructed in analogy to Ref. ^^ilSj and combined 
with the contour deformation prescription of Ref. [7]. In addition, a new treatment for the 
UV divergences is introduced. As a ffist step, it is shown that this approach is already 
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Figure 1: Topology of a general one-loop integral. 



useful at the one-loop level, leading to numerical integrals with good convergence behavior. 
Secondly, it is demonstrated how it can be applied to two-loop diagrams. In the present 
paper, only two-loop integrals with IR singularities in one of the subloops are considered. The 
extension to overlapping IR singularities in both subloops will be delegated to a forthcoming 
publication [H], but it is relatively straightforward since most of the required subtraction 
terms can be constructed from the results of Ref. jl5j . 

A public computer code for the evaluation of arbitrary one-loop integrals is provided, 
based on the method presented here. It is planned to incorporate two-loop cases into future 
versions of this code. 

In the next section, the method for the numerical evaluation of one-loop integrals, using 
subtraction terms and contour deformation, is explained. It is demonstrated how it can 
be cast into an algorithmic form and its application in several examples is demonstrated. 
In section [3] the extension to the two-loop level is discussed and illustrated through a few 
examples. Finally, the main results are summarized in section |H 



2 One-loop integrals 

A general one-loop integral, see Fig. [1], may be written as 

D^'\k) = [k' - ml] [{k - pif -ml]--- [{k - p^f - ml], (2) 

where N{k) is a polynomial in the loop momentum k. This integral may have soft, coUinear 
and ultraviolet divergences. The first two are canceled in the integrand by using suitable 
subtraction terms, similar to Ref. [12] (see also Ref. [I3]). However, no UV subtraction terms 
are used for the one case here; instead the UV poles are evaluated explicitly after introducing 
Feynman parameters. 
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2.1 Infrared subtraction terms 

ropagat( 

{pi~i —piY = ^l-ii N{k) 7^ for k — > pi. The soft subtraction term is given by 



The zth propagator of integral ^ produces a soft singularity if rrii = 0, {pi+i —piY = rn'^+i, 
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where, for later convenience, the last line has been written in terms of the common denom- 
inator of eq. ([2]). The integrated soft subtraction term reads 

1 
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where e = {4 — d)/2, s = {pi+i — Pi-iY + ie, and 



j— 1, j, i+1 



{pi -PjY - 



A coUinear singularity singularity is encountered if = mj_i = 0, {pi — Pi-i)^ = 0, and 
N{k) 7^ for k — )■ pi. The simplest subtraction term for this singularity is given by 



^(1) 

{k - p,^,nk - p,y ' ■ 11 ^ - p,y - 



N{k=p,) H 



(A; - - 

{Pi-PjY -m] 



(9) 
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In dimensional regularization, the integrated coUinear subtraction term is simply zero, 



I dk G« = 0. (10) 

In other words, by subtracting G^oii ^oi^ the integrand the collinear singularity is trans- 
formed into a UV singularity, which will be handled as explained below. Note that this 
collinear subtraction term only works when applied to physical amplitude, i. e. to a gauge- 
invariant set of diagrams [T2] . 



2.2 Variable mapping and contour deformation 

After subtracting the IR singularities as described above, one arrives at a loop integral that 
only contains UV divergences, 

- - E/ dkGi^, TkGil, (11) 

where the sums are over all soft and collinear singularities in It has a structure similar 
to eq. ([2D: 

Til) - fit ^'-^g(^) (^O) 

By introducing Feynman parameters and shifting the loop momentum, this expression can 
be cast into the form 

/« = jdx,... dx^ 6{1 - Ex.) I dk j^^, (13) 

where N{k) is a polynomial in k and in the Feynman parameters xi,...,Xn, while A is a 
polynomial in xi, ...,Xn- It is convenient to map the Feynman parameters onto a hypercube: 

xi = l- 1/1, 

X2 = yiil-y2), I dXi...dXnS{l-Y.Xi) 

\ ^ .1 (14) 

Xn = yi- ■■yn-2yn-i, 

The /c-dependence in the numerator of (fT3ll is eliminated by using the tensor reduction 
formula 

dk 



[A;2 - A]^ 

1 ^ p (15) 

r\\d{d + 2)---{d + r-2) ^ ^ ' J [k^ - A]-' 



permut. 



where the sum is over all permutations of /ii, Thus one arrives at 
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dyi... dy, 



n-l 
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where the Cj are polynomials in ?/i,...,?/„ but independent of k. Integrating over k, one 
obtains 



/(I) 
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r(n- 1) 

which can be expanded in e, so that the UV singularities appear as 1/e poles. Thus the 
integral takes the form 
»i 



reg 



dyi... dyn-1 



Do + log(A - «e)) +Di{A- te)-^ + Di^A - te] 
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1^ 



where the ie from the propagators have been made explicit again. A and the Di are polyno- 
mials in the variables yi, ...,?/„_i. 

The coefficients in the e-expansion in eq. f|T8l) are finite, so that the integration over 
the parameters j/i, can be performed numerically. However, the polynomial A{y) 
can have zeros inside the integration region, which happens if the loop integral has physical 
thresholds. While these singularities are formally integrable, they lead to problems for the 
numerical integration. A solution is the deformation of the integration contours into the 
complex plane, by using the variable transformation 



yi = Zi - lAZiX i - Zi) — , Q < Zi< 1. 
To leading order in A this produces a negative imaginary part in A: 

A{y) = A{z) - a 5^ z.(l - z,) + 0{\'). 



So if A is chosen small enough, the integral 



reg 



dz\ . . . dz^_i 



d{yi,...,yn-i) 



d{zi, . . . , Zn-l) 



Do(^- + log{A-ie)^+DiA-'+D2A-^ + . 



(19) 



(20) 



(21) 



is well-behaved for numerical integration. In most cases, a good choice for A is roughly 0.5. 



2.3 Numerical examples 

To demonstrate the fiexibility and efficiency of the method described in the previous sec- 
tion, several example calculations have been carried out and, where applicable, compared to 
existing analytical results. 

The numerical integration has been performed with the VEGAS and Cuhre algorithms 
of the Cuba 1.4 library [IB]. Timing information is given for running on a single core of a 
Intel® Xeon® X5570 processor with 2.93 GHz. 
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Figure 2: Results for the four-photon amphtude A1(77 — )■ 77) (normahzed relative to 
o?) for different photon helicity combinations. The dots correspond to the results obtained 
with the numerical method presented here, while the lines depict the analytical result from 
Ref. The contour deformation parameter has been chosen A = 1. 



Figure 3: Feynman diagrams for the one-loop QED corrections to Ueiy^ e /i^. 

1. Let us start with the amphtude for two-photon scattering 77 — )■ 77 at the one- loop level, 
with electrons running in the loop. This process does not involve any physical singularities 
and thus is a simple test of the contour deformation and numerical integration. The six 
contributing diagrams are grouped into three groups, according to the structure of the loop 
denominators. Each group is evaluated separately, using the Vegas algorithm with 10^ 
integration points each, which in total takes 4.8 s evaluation time. Results for different 
helicity combinations of the photons are depicted in Fig. [21 and compared to analytical 
results from Ref. [21]. The numerical integration error ranges between 0.2% and 0.5%; error 
bars are included but not visible in the plot. 

2. The one-loop QED corrections to I'e'^fj. e~^~^ is an example with UV, soft and collinear 
singularities. The six diagrams in Fig. [3] are combined into a single expression with one com- 
mon denominator, from which the soft and collinear divergences are subtracted as explained 
in section [2TT1 
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this work 


analytical 




-0.625 


-0.625 




1.13113365 - 1.96349541i 


1.13113365- 1.96349541i 


0(6°) 


3.27791(4) + 1.83014(4)i 


3.27792343 + 1.83014477i 



Table 1: Numerical results for the next-to-leading order QED corrections to the squared 
matrix element of Ve^^i e~n~^. Specifically, the table shows numbers for 2 Re{Ai^'^^ Ai^^^} , 
for the input parameters Mw = 1, s = 2, t = —1, and A = 0.5. The results from the 
numerical method described in this paper are given in the second column, together with the 
numerical error from 5 x 10^ integration points of the CUHRE integration routine. For com- 
parison the third column lists the results obtained with the traditional analytical approach 
using Passarino-Veltman reduction. 



The results presented in Tab [T] correspond to the next-to-leading order correction to the 
squared matrix element, given by 2Re{A^(°)A^(^)}, where A^(°) and M^""^ refer to the tree- 
level and one-loop QED corrected matrix elements, respectively. The numerical integration 
has been performed with the Cuhre algorithm, using 5 x 10^ integration points (resulting in 
2.4 s running time). Note that the e"^ pole stems solely from the soft subtraction term and 
thus is known analytically. For comparison, the table also shows the result from a traditional 
analytical calculation using Passarino-Veltman reduction [12] and explicit expression for the 
well-known basic integrals, see e.g. Ref. [20] . 



3. A third example, with a more difficult denominator structure, is the scalar hexagon 
integral depicted in Fig. |H It contains several mass scales, has three soft singularities, and 
requires contour deformation. Numerical results are shown for the input values 



(22) 



ml = 1.0, 


Pi = 


(0,0,3), 


P2 = 


(0,0,-3), 


ml = 0.25, 


P3 = 


(0.5,0,0.5), 


Pa = 


(-0.2,0,0.1), 


ml = 4.0, 


Pb = 


(0,1.37626,-0.3), 


P6 = 


(-0.3,-1.37626,-0.3), 



and A = 0.5. Using 10'' integration points with the Cuhre algorithm one obtains 

0{6-^) : - 0.0044718804 + 0.0120697975i, 
C(£°) : - 0.04383(14) + 0.00790(14)?, 

which takes 45 min to evalute on the test computer used here. 



(23) 



2.4 The computer package NICODEMOS 

The procedure outlined in sections 12.11 and 12.21 has been implemented into the public com- 
puter code NICODEMOS (Numerical Integration with COntour DEformation and MOdular 
Subtractions). The package contains a Mathematica module, which performs the appli- 
cation of subtraction terms and Feynman parametrization (with or without deformation). 
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Figure 4: Scalar hexagon diagram with soft singularities, labeling the masses of all internal 
and external lines and the external momenta. 

The user needs to supply the input expression and the location of IR singularities. The 
Mathematica code then produces a Fortran executable, which performs the numerical 
integration for a given set of numerical input parameters. 

NICODEMOS is available for download at http://www.pitt . edu/~af reitas/ and can be 
used freely, provided the source is acknowledged and properly cited. Version 1.0 can only 
handle one-loop integrals, but it is planned to expand the functionality to include two-loop 
integrals in future versions. 




3 Two-loop integrals 

In this section the extension to two-loop diagrams is discussed. The two-loop integrals may 
contain UV singularities as well as IR singularities in one subloop, while the case with IR 
singularities in both subloops will be addressed in a future publication [H] . 

Of particular interest are relatively complex two- loop diagrams, i. e. diagrams with a 
relatively large number of external legs. Therefore the following discussion will not consider 
some special cases that occur only for two-loop tadpole and selfenergy diagrams, since these 
can be evaluated with existing methods, see e.g. Ref. [22] • 

Following the notation of the previous section, a two-loop integral is given by 

D^'\h, h) = [kl - ml][ih - Pif -m\]--- [{k, - prf - ml] (25) 
X [{k2 - Pr+if - ml^^] ■ ■ ■ [(^2 - Psf - ml] 
X [{ki -k2- Ps+i^ - "^s+i] ■ • ■ [{ki - /C2 - Pn^ - ml], 

where N{ki, ^2) is polynomial in ki and k2 and in the external momenta. 

If one of the two subloops has an IR singularity it can be subtracted analogously to the 
one-loop case. For example, the scalar diagram in Fig. |5] 



^fitE]= / dkidk2 (26) 

1 



X 



[ki — m^] [(ki — pi)^] [(/ci — P\— P2)^ — m^] [(A;2 — Pi)^ — m^] \{kx — A;2)^ — m^] 
has a soft singularity for k^ — > Pi. The subtraction term is constructed as in eq. (|3]), viz. 
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Figure 5: Scalar two-loop diagram with soft singularity. The Ggures also labels the masses 
of all internal and external lines and the external momenta. 



This expression factorizes into two one-loop integrals, which can be integrated analytically. 

In the same way, one can subtract collinear divergences in one subloop with the subtrac- 
tion terms introduced in section 12.11 



3.1 Ultraviolet divergences 

Since two-loop integrals may have overlapping UV singularities from both subloops, the 
corresponding 1/e poles cannot be computed directly as in eq. ( !T8|) . Instead one has to 
introduce subtraction terms also for the UV divergences. The UV subtraction is performed 
in two steps: 



1. The global UV singularities of both subloops can be obtained by performing a Taylor 
expansion of the two-loop amplitude in terms of the external momenta. For all physical am- 
plitudes with three or more external legs only the leading term in this expansion contributes 
to the global UV divergence. Two-loop tadpoles and selfenergies are not considered here, as 
mentioned above. Thus the global UV subtraction term is defined as 



^(2) 
*-^glob 



N{h,k2 



(28) 



P,=0 



The integrated subtraction term J dkidk2 G'^^I^ consists of two-loop vacuum integrals, which 
can be evaluated analytically with the methods of Ref. [23] . 



2. The remainder Igs = iLg — J dkidk2 G^^^^ can still contain a UV singularity in one of 
the subloops, or in both. The latter case, however, only occurs for tadpole and selfenergy 
diagrams and thus will not be considered here, since it can be handled more efficiently with 
other methods [22] . 

Let us then assume that only the subloop with loop momentum ki has a UV divergence. 
The ki integral is now turned into a Feynman-parameter integral by following the steps in 
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Figure 6: Typical diagram for the two-loop correction to the Zbb vertex. 



eqs. ffT51) - ffTB]) . leading to 



r(2) 

gs 



dyi... dy, 



dkidk2 



Co 



[kl-A] 



m—l 



(29) 



where the Ci and A depend on k2 and the Feynman parameters yi, ym-i- Here m = n+r—s 
is the number of propagators with ki. Note that the Ci are rational functions containing the 
propagators that depend on ^2 only. At this point, a sub loop UV singularity is indicated by 
a term of the form Cj/[kf — A]^ in eq. ( !29|) . All higher powers of [kf — A] in the denominator 
are UV-finite. The subloop UV divergence can be subtracted by 



G. 



(2) 
sub 



dyi... dy, 



m— 1 



[kl - /. 



212 ■ 



(30) 



where /i is a suitably chosen mass parameter. The integrated subtraction term factorizes 
into two one-loop integrals: 



dkidk2 G 



(2) 
sub 



2-e 



dyi...dym-i / dk2Cj. 



(31) 



The one-loop integral over k2 can now be carried out with the procedure of section 12.21 

The remaining two-loop integral /rem = /reg — / dkidk2 Cgj^j^ — / dkidk2 Gg^^ is finite. Af- 
ter introducing FejTiman parameters for the ^2 integral, which are mapped onto a hypercube, 
and shifting the /c2 loop momentum one arrives at an expression of the form 



j(2) 



dyi... dyn-i / dkidk2 



Cl{k2 



+ 



C2{k2 



[akf + Pkl - AY [akl + Pkl - A] 



n-l 



+ 



(32) 



where a, (3 and A are polynomials in the Feynman parameters yi, ...,yn-i, and Ci{k2) are 
polynomials in k2 and the Feynman parameters. Now one can perform the tensor reduction 
for ^2, the ^2 loop integration, and — if necessary — the contour deformation in analogy to 
eqs. (115!), (inD, and 



3.2 Numerical examples 

1. A typical example with global and subloop UV singularities is given by the diagram in 
Fig. [6], which contributes to the two-loop corrections to the effective weak mixing angle of 
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this work 


Ref. [2] 




-2.30183413 


-2.30183413 




5.07108758 


5.07108758 




8.326(1) 


8.3259 



Table 2: Numerical results for the contribution of the diagram in Fig. to the effective 
weak mixing angle sin^^gg, in units of (a/47r)^. The input parameters are chosen as in 
Ref 0/; Mz = 1, Mw = 80/91, mt = 180/91. The deformation parameter is set to X = 0. 
The numbers in the second column have been obtained using 10^ integration points of the 
CUHRE integration routine, with the integration error given in brackets. The last column 
shows the results from to Tab. 1 in Ref. ^2] for comparison. 





this work 


Ref. [15] 




-0.43040894 + 1.40496295i 


-0.43040894 + 1.40496295i 




-3.53105702 


-3.53105702 


0{e') 


-1.93471(1) - 2. 08763(l)i 


-1.93471213 - 2.08762578i 



Table 3: Results for the numerical evaluation of the scalar diagram in Fig. O The second 
column shows the values obtained with the method of this paper, using m = l,s = 5, A = l 
and 10^ integration points of the CUHRE algorithm. For comparison, the last column gives 
the analytical result from eq. (343) in Ref. ji5| . 

bottom quarks, sin^ 9^^. Since this diagram does not have any physical cuts, no deformation 
of the integration contour is needed and one can choose A = 0. Tab. [2] shows the numerical 
results obtained with the algorithm of the paper, using 10^ integration points of the Cuhre 
routine. The final answer is the sum of the global UV subtraction term, the subloop UV sub- 
traction term and the remaining finite two-loop integral. The integration error is negligible 
except for the finite 0{e^) term. The numerical integration takes about 13 min to evaluate 
on one core of a Intel® Xeon® X5570 processor with 2.93 GHz. For comparison, the 
table also shows the result of Ref. [2] , which has been obtained with the Bernstein- Tkachov 
(BT) method. The two results agree very well within numerical integration errors. 

2. Fig. [5] shows an example of a scalar diagram with both UV and IR divergences. The soft 
subtraction in eq. (1271) is used, but its integrated form needs to be evaluated to order 0{e), 
due to the presence of the UV singularity. Numerical results obtained with this method are 
shown in the second column of Tab. [3l using 10^ integration points of the Cuhre routine, 
corresponding to a running time of 2.6 s. For this particular diagram an analytical result 
has been obtained previously in terms of harmonic polylogarithms [15j, which is also shown 
in the table for comparison. 



11 



4 Summary 



This paper presents a procedure for numerically computing one- and two-loop integrals using 
subtraction terms for the singular pieces of the integrand. A set of subtraction terms for the 
removal of ultraviolet, soft and collinear singularities has been described. The subtraction 
terms themselves can be easily integrated analytically in dimensional regularization, while 
the difference between the full loop amplitude and the subtracted contributions is finite 
and can be evaluated numerically after taking the limit to four dimensions. The approach 
is apphed to general one-loop cases, as well as two-loop integrals with global and subloop 
ultraviolet divergences but infrared divergences in only one of the two subloops. The exten- 
sion to overlapping infrared singularities in both subloops will be discussed in a subsequent 
publication. 

The finite numerical integral (after application of the subtraction terms and expansion in 
the integration dimension) is evaluated in Feynman parameter space. The integrand can still 
have singular points inside the integration regions if the diagram has physical thresholds. 
These singularities are formally integrable but are problematic for numerical integration 
routines, so that they must be avoided by deforming the integration contour into the complex 
plane. 

The usefulness of the proposed procedure has been demonstrated with several one- and 
two-loop examples. One obtains an overall good convergence behavior of the numerical inte- 
gration, although problems can occur for integrals with pinch singularities, which typically 
correspond to threshold configurations. It may be possible to regularize the pinch singular- 
ities with additional finite subtraction terms, but a detailed investigation of this matter is 
left for future work. 

The algorithm presented in this paper has been implemented in the public computer pro- 
gram NICODEMOS (available at |http: //www.pitt . edu/~af reitas/ ). It is based on Math- 
EMATICA for symbolic manipulations and produces a Fortran executable for the numerical 
evaluation. The current version 1.0 only contains one- loop functionality, but an extension 
to the two-loop level is planned for future releases. 
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